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ABSTRACT 

A brief review is presented of various problems which are con- 
fronted in the development of an unsteady finite difference poten- 
tial code. This review is conducted mainly in the context of what 
is done for a typical small disturbance and full potential method. 
The Issues discussed Include choice of equations, linearization a nd 
conservation, differencing schemes, and algorithm development. A 
number of applications, including unsteady three-dimensional rotor 
calculations, are demonstrated. 
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ONE OF THE MOST IMPORTANT REALIZATIONS in fluid mechanics research was 
that transonic flow, for all its apparent complexity, is largely 
describable by potential theory. This fact was clouded by a thicket of 
problems concerning tunnel turbulence, wall and scaling effects, and sepa- 
ration. We now know that many of these problems are magnified by the 
inherent susceptibility of the inviscid transonic flow to unsteadiness. 

Of course, basic flow researchers were undoubtedly prejudiced in favor of 
potential theory from the outset because the inviscid, irrotational 
approximation is a tremendous simplification. Nevertheless, real progress 
did not occur until the explosive development of finite difference methods 
which occurred in the 1970 *s. 

Spurred on by the promise of more cruise-efficient transport air- 
our greatest progress has been in steady flow prediction. However, 
unsteady transonic flows are important because the aircraft structural 
response can induce large and heretofore unpredictable shock excursions. 

The flow on helicopter rotor tips is an even more interesting example of 
unsteady transonic flow. In this case, unsteadiness is not only forced by 
structural deformation but also by a free-stream flow which rapidly varies 
in speed and direction and contains large wake disturbances from previous 
blades. And so, unsteady transonic flow remains a rich mine of Important 
problems waiting to be solved. Potential methods will probably play a 
dominant role in this process — not only because of validity but also 
because it promises to be the most efficient method. 

Efficiency is a more Important matter for unsteady computations than 
for the steady case. This is because we require the resolution of physical 
time and do not have the benefit of acceleration methods that are used in 
steady problems. In general, however, the unsteady and steady methods are 
much the same. At present, the most versatile and efficient unsteady 
methods are two- and three-dimensional small disturbance codes, but full 
potential methods are currently under rapid development. The following 
discussion will review many of the Important algorithm and code develop- 
ment issues in the context of small disturbance and full potential methods. 

FORMULATIONS OF THE PROBLEM 

The starting point for the various potential formulations is the 
mass conservation and Bernoulli’s equation (shown here for an inertial 
reference frame) 


+ V • (p7(J)) = 0 

1 

(1) 


(2) 


These equations have the advantage (when combined) that only one flow vari- 
able, <P, need be solved for. However, these equations are only an approxi- 
mation to the exact inviscid equations because they assert that mass, 
energy and entropy are conserved throughout the flow field. The component 
of momentum normal to any shock is not conserved in this approximation 
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(1,2)*. This is a notable difference from the Rankine-Hugoniot equations 
where shock drag comes about through the change in entropy across the 
shock. The error thus induced is generally not excessive until shock 
strengths are attained that involve separation and the necessity to abandon 
the inviscid approximation. Conservative formulations are required to 
closely approximate Rankine-Hugoniot results. Nevertheless, nonconserva- 
tive formulations are still actively employed, due to the happy accident 
that errors of conservation frequently have an effect similar to boundary- 
layer corrections on shock location. 

In attacking a given problem it is first necessary to express Eqs. (1) 
and (2) in the relevant body-fixed coordinate system for the problem to be 
solved. The transformation (3) which encompasses both wings or rotors in 
edgewise motion is 


r' = U t + R(t)r 

00 

t' = t 


(3) 


where r' = (x’,y',z') and r = (x,y,z) are the Inertial and body-fixed 
coordinates (see Fig. 1), respectively, and R(t) is the rotation matrix 

cos SI t -sin fit 0 

R(t) = Isin H t cos t 0 

L 0 0 1 

(Note that for no rotation, = 0, Eq. (3) reduces to the usual Galilean 
transformation.) Under this transfoinnation, Eqs. (1) and (2) become 


Pj. + V • jpvj (4) 



where a = U^(i cosHt - j sln^t) - x r is the undisturbed free-stream 
velocity seen by an observer in body-fixed coordinates and 


V = a + 7(|. (6) 

is the flow velocity seen in the moving coordinates. A common misunder- 
standing concerns the validity of potential theory wljere rotary motion is 
involved. In fact, the motion of the coordinates, -a, has no bearing on 
whether or not a potential exists. However, since a need not be an 
Irrotatlonal function, V is not generally expressible as the gradient of 
a full potential. Rather, defines a disturbance about a, as seen in 
Eq. (6). 


* 


Numbers in parentheses designate References at end of paper. 
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Equations (A) and (5) are here written for a generalized, moving coor- 
dinate system C»n»?,T where 


giving 

where 




n 

? 

T 

3^(p/J) + 3^(pU/J) 
U - + VC • 

V = + Vn • 

w = + 7? • 


“ C (x',y',z',t) 

“ n (x',y',z',t) 

= ? (x',y',z',t) 

=> t 

+ 3^(pV/J) + 3^(pW/J) 
((J»^VC + 

(<(i^7C + <J»^7q + <j>^7c) 
(4,^7C + <J»^7q + <fr^75) 


= 0 


(Aa) 


are here the contravariant velocity components and have the form of Eq. (6), 
J is the Jacobian 1 3(C,q,c) /3(x* ,y' ,z ') | and 7 is the cartesian 
gradient operator. In the same manner, Bernoulli's equation becomes in 
these coordinates 


p = {l - 

^ (5a) 

+ (C (f_ + q <J) +5 + ?<f> )^l ^ ^ 

y C y q y^c '“^Z^C 'z^q ^Z^C jf 

where we use the nondimensionalizations p p/p„, x x/i, y -*■ y/£, 

^ ^ z/£, $ ->■ 4)/£c^, and x tc„/ly and the tilda is suppressed. Here, £ is 
a reference length such as the airfoil chord. 

SMALL DISTURBANCE EQUATIONS— The classical small disturbance deriva- 
tion involves substituting Eq. (5) into Eq. (A) to eliminate p. The 
resulting equation is nondimensionalized, scaled, and higher order terms 
are eliminated, subject to the limit process that (1 -m 2)/62/3 = 0(1) as 
M •+■ 1 and 6 -► 0. The resulting equation takes the form, 

Aij*.... + Bij> ^ = F + <ii + c<i) + Dd) ( 7 ) 

tt ^^xt X ^zz ^yy ^xy ' ^ 

where 

^ (J>/U^£62/3 

0 = frequency of unsteady motion (the rota- 
tion rate for a rotor) 

M = a characteristic Mach number 

U^ = characteristic speed (V^ for a wing, 

$1R for a rotor) " 

AR = R/£, aspect ratio 
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k = fi£/Uc, reduced frequency (for a rotor 
k = 1/AR) 

6 = blade thickness ratio 
A = M2k2/62/3 
B = ZM^kf 762/3 
C = 1/AR262/3 
D = Bg 

f = y + pcos t (for rotor), or 1 (for a wing) 

g = X + sin t (rotor) , or 1 (wing) 
y = J2R/V , rotor advance ratio 
t « nt' 

X = xWi 

y « y*/R (R is either a rotor radius or the 
wing span of a fixed wing) 


I = chord 

Equation (7) is not unique and an assortment of modifications and 
additions have been made to improve its ability to handle oblique shocks 
(see Ref, 4). Nevertheless, they all have the same general form and have 
essentially identical unsteady terms. 

The pressure coefficient in the small disturbance approximation is 
given by 

P - P 

c = J 1 - - 2f62/3 U + k6 ) (8) 

ipu/ >= ^ 

The body surface boundary condition Is the familiar slope condition 
applied at a mean surface. 


'^z “ ^^^^t V 

where the body surface Is described by z = 6b(x,y,t). 

There Is an additional boundary condition that pressure be continuous 
across the stagnation streamline behind and on the trailing edge of an 
airfoil (Kutta condition). Using Eq. (8), this condition Is expressed as 
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° ^ = [[^]] ( 10 ) 

This equation expresses the downstream convection of vorticity (or the 
potential discontinuity). It happens that (p Is also discontinuous in 
the airfoil wake. This is seen by applying Eq. (7) (in linearized form) 
across the wake and substituting Eq. (10) to obtain 

X ■ ■ - [“ * - f)] ■'zx - "yy <11) 

where 


This quantity, x» is zero only in steady two-dimensional flow or steady 
nonlifting three-dimensional flow. 

Finally, finite difference methods require far-field boundary condi- 
tions. Typically, one specifies some combination of ({> = 0, (ji = 0, or 
Cp = 0 on the outer boundaries and subsequently relies upon enlarge 
boundary distance (often over lOOi) to dissipate the ensuing wave reflec- 
tions. However, it has been shown that good approximate nonreflecting 
boundary conditions can be constructed (5—8) . For example, consider 
Eq* (7) in its two-dimensional, linearized form 

The wave information for a plane wave, $ - glCwt+5x+tz) ^ satisfy 

Aii)2 + Bo)C = C52 + ^2 


A condition which prevents reflection of upstream-moving waves at the 
front grid boundary is 


^ „ -Bg + ✓(B^-f4Ag)£^ - 4^ 


This expression does not transform back to a simple local differential 
operator. However, in the limit n/g 0 (for wavefronts parallel to the 
upstream boundary) we obtain the expression 


lx) 


( 


✓B2-f4Ao - B 
2A 


■) 




whose inverse Fourier transformation yields 


. „ /B2-KAg - B . 

“^t 2A ***x 


Similar expressions can be obtained for all boundary faces. 


( 12 ) 
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Although its limitations are many and well known, the small distur- 
bance approximation (with judicious use) often gives excellent results and 
remains in active use. Perhaps the greatest importance of small distur- 
bance theory is that it contains (in simplified form) all the important 
issues contained in the full equations and hence constitutes a good test- 
ing ground for future techniques. There is, perhaps, one exception to this 
statement. It concerns the fact that when Eqs. (4) and (5) are combined to 
produce a single equation for <j>, a nonconservative formulation results. 
This problem does not occur in classical small disturbance theory. 

FULL POTENTIAL EQUATIONS~For computational efficiency, we usually 
seek to formulate the potential equations so that we need only one depen- 
dent variable. This is easily done, since Eq. (4) can be rewritten as 


d2_ 

3$ 


where 


3'^ > f . 3^^ . 3^A . . 9p 3$ 
3t + ST 


+ $ 


9p 9^ 

94^ 9y ^ 


I. _f>2-Y 


3$ 




3t 


_3 

3x 




3y 


_3 

3z 


$ ^ 


3$ 


li 

3z 


) 


(13) 


is a noncoramuting differential operator. For this equation we assume the 
nonlinearizations pertaining to Eq. (5) (for the previous and following 
equations we assume pure translational motion) . On applying the isen- 
troplc relations we obtain the familiar textbook form of the potential 
equation 


^t 


2u^ ^ + 
xt 


2v$ 


yt 


a (C^-U^)$ + - 2$ $ $ 

XX yy 2Z X y xy 

- 2 $ $ $ — 2 ^ $ $ ( 1 ^) 
X z xz y z yz 


The problem arises, however, that equation (14) cannot be brought back 
into conservation law form with (j> retained as the dependent variable. 

In examining Eq. (4) we see that p presents no problem in the 
special derivative terms because it can be evaluated from the previous 
time step (that is, it is taken as the leading term in a Taylor series in 
time) . The required implicit spacial tertos involving ^ then come about 
through V. What is required is a way to obtain a conservative temporal 
function of $ from the p term. This can be accomplished through the 
following linearization, 

where the subscript, o, denotes a nearby known state or solution. With 
this expansion the density time derivative becomes 

"t * ^ ((100*) * A (‘’o-(H)„«o) 

which provides a conservative time-differenced function of $ (9,10). 


( 16 ) 
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The above linearization can be avoided for a conservative formulation 
only if one is willing to retain two dependent variables. This approach 
huM h.M.n l.'ikcn by ailpniai. aa.l .lameson (II) wl.o .solve riie nrst-orllcr 
MyMi.sM roiiMlfii Ing o| Uir ma.MM c.imorvaL Ion f(|unLloii and IUtiiouJU'h oaua- 
tion rewritten in the form 


h(p,74>,c^) = 0 


SPACIAL DIFFERENCING 


The primary stability issue in transonic finite difference problems 
(steady or unsteady) concerns proper spacial differencing in the subsonic 
and supersonic flow regions. Consider first the simple steady two- 
dimensional model problem 


(1-MJ) + 4>yy = 0 (17) 

and let the difference operators , and A be defined as 

- (4>i + (p^^i)/^x, A^(l> « (<l>i+i - (t>±)/AXy etc. ^It is well known that the 
following difference schemes 

(1-mJ) V^A^4> + 7yAy(fr = 0 , < 1 (18a) 

Vx* ° ^ ^ ( 18 b) 

are, respectively, suitable for the above model elliptic and hyperbolic 
problems. The differencing equation (Eq. 18b) is divergent if < 1 
while (Eq. 18a) is convergent for M« > 1 only if | Ax/[ (l-Mo?) Ay] I < 1, 
an impractical restriction as Mob -»• 1. Murman and Cole (12), aware oT 
these constraints, introduced type-dependent differencing for the transonic 
small disturbance equation. In this approach the streamwlse spatial dif- 
ference operator is either central or backward, depending on the local 
Mach number. 

In order to obtain a stable conservative streamwlse difference opera- 
tor for Eq. (7), it is first necessary to time-linearize the nonlinear 
flux term F. One way to do so (13) is 


j.n+1 0= f” + 



This time-linearized flux is differenced and defined at each midcell of 
the computational grid as 



n 

i+h 


+ 



+ 



7j.(,|>"+l-<^”) 
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wnere / n 

= a - m 2 |(Y+l)fA^(j)" + k(Y-l)7^())“J 

and 

(^)” ■ 

For this flux operator Murman's conservative switching scheme (14) is 

* V = S ‘^l-%-’'i-3/2)] (15) 


[0 , < 0 

= a - m2(y+D<|)^ - M2(Y-l)k(j)j. 

Note that for greater time accuracy, the Crank-Nicolson time averaging of 
spacial operators^^l often applied. For instance, we could define the 
flux term F = (F + F^)/2 and subsequently apply the above lineariza- 
tion and switching scheme. 

The nonconservative full potential equation (Eq. 14) can be type- 
differenced in a manner very similar to that of the small disturbance 
equation. Equation (14) can be rearranged in the cannonical form (15,16) 


where 


and 


♦tt 5- (!' 

where s is a coordinate locally aligned to the stream direction and n 
is a normal coordinate. For the two-dimensional case 


*^nn 


(v2(|) +2uv(}> +u2(|) ) 

q2 ^xx ^xy ^yy' 


and 


|> “ (u2(|. +2 4) +v2(j) ) 

'^ss „2 ’^xx uv’^xy ^yy 


2q<|.gt = 2u<fr^t 


It is seen that the steady part of Eq. (20) has the same form as Eq. (18). 
Therefore a stable scheme results when (fun 4® always central-differenced 
and <fss 1® type-differenced based on the sign of c^ - q2. 

This assertion of stability, by analogy to known stable methods, can 
be applied to the full potential equation. For this case we require the 
following difference scheme for the model equation, Eq. (17). 
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This scheme, though less accurate than Eq. (18a), is stable for the 
entire Mach number range with no need to switch difference operators. 

Now consider the spatial derivative term axPifix of Eq* (4). A 
local linearization of this term gives 

[(‘’♦x)„ + (♦xlt + P 

which upon substitution of Eq. (13) and assumption of steady small distur- 
bance flow yields the approximate expression, 

^x^P'*'x^ “ ^x^Po*^x^ - (23) 

If we further assume and (u/c)^ to be spacially constant we have 

^x'»*x> - 0 [♦xx-“''+xx] (24) 

which demonstrates an approximate equivalence between the full potential 
equation and the model Eq. (17). The Importance of Eq. (24) is that the 
origin of the second on the right-hand side is in the evaluation of 

p from Bernoulli’s equation. Therefore, the differencing scheme of 
Eq. (21) is significant because it is roughly equivalent to evaluating 
^xP'J’x with a centered scheme employing an upstream biased p. The sta- 
bility of this density biasing has been demonstrated in Refs. (17-19). 

This density biasing will be expanded upon in the following section. 

NUMERICAL ALGORITHMS 


The final step in a finite difference solution scheme is the effi- 
cient solution of the system of algebraic equations which result from 
the various discretizations. Because these systems are all too large to 
be solved efficiently (in spite of their sparseness) , it is necessary to 
reduce them to a more manageable form. All the unsteady schemes today 
use some sort of approximate factorization. That is, the system matrix 
is replaced by a product of easily solved submatrices. In general, this 
product is not equal to the original system matrix. However, it is pos- 
sible to keep the discrepancies within the bounds of the discretization 
error . 

SMALL DISTURBANCE EQUATION-The germ of the approximate factorization 
idea is quite old, having been first expressed in the ADI method. A 
typical ADI scheme for Eq. (7) is 


Step 1. V (^-(j)”) = D F(^) + C7 A (()” + 7 A <|>” + D7 

4i*-A X yy zz y 


Iv" 


D < 0 

b > 0 
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step 2. = D^F($) + + D7^ 

Step 3. A7^7^^"'^^ + 37^7^4."'^^ = D^F($) + CF^Ayl + 


7y^ , D < 0 
Ay^ , D > 0 ’ 


+ D7 ' 

X 


Vyl 

Ayl 


D < 0 
D > 0 


Note that the cross derivative term is treated implicitly and is always 
an upwind difference (20,21). The advantage of this ADI scheme is that 
each step involves a simple matrix inversion (steps 2 and 3 require a 
tridiagonal and step 1 requires a quadradiagonal inversion) . Note also 
that these steps are written so as to maximize their similarity to Eq. (7). 
However, they are never actually solved in this form. It is more effici- 
ent to subtract steps 1 and 2 from steps 2 and 3, respectively, and solve 
the resulting equations. The insertion of the <|)tt discretization in 
the above algorithm is very convenient and was first introduced by 
Rizetta and Chen (22) . 

It can be shown that the above ADI scheme constitutes a product of 
temns which are a good approximation to the original difference equation. 
This scheme was originally devised with the idea that each step should be 
a consistent approximation of the original difference equation. However, 
if this constraint is relaxed there are a number of factorizations which 
can be found. To see how this is done, consider the model equation 


^xt ^xx ^yy 


which is differenced as 


V 7 = 07 A + 7 A (for 0 > 0) 

t X X x'^ y y 


and then rearranged as 


(7 -07 A -7 A = 7 <j)" 

X X X y y X 


The left-hand side of this equation can be approximated by the product 
(l-0Ajj) (7jj-7yAy) if the error term 07xAy7y is eliminated. This elimi- 
nation is easily accomplished using the known information at step n. 

And so, an appropriate approximate factorization (commonly called the AF2 
scheme) is 


(l-07^)(7^-7yAy)<|>“'^^ - (7^+0A^7yAy)4." 

This gives rise to a two-step procedure — one step being a bi-diagonal 
and the other a tri-diagonal inversion. A further simplification is to 
solve for ((J»’^"*‘l - (J)!') rather than (j)n+l in order to simplify the evalua- 
tion of the right-hand side. 



11 


CONSERVATIVE FULL POTENTIAL METHOD-The above factorization approach 
will now be applied to the unsteady full potential equation. We will 

® simplified form of the equation for two dimensions. Con- 
sider Eqs. (4) and (5) in their two-dimensional form subject to the 
^pping C =» C(x),n = n(x) and assume a steady free-stream motion allow- 
ing (p to be expressed as <|> -»■ M^x + (j) . This gives 

\(P/J) + ♦ J = 0 (25a) 

P “ 1 1 - [♦x )]) <25b) 

On application of the linearization of Eq. (15), Eq. (25a) is differenced 
as 






(26) 


where 3 = P^“Y, implies division by J, and the reference state 0 is 
now taken to be the previous time step n. The special operators in 
Eq. (26) are defined as 


^ ^i+1 ~ ^ ^i-i 


6 

n 


( >!+!-( >i-l 


(27a) 

(27b) 





[' 




Pj+1+ Pj 
2 


+ V 


i+1 


(l+e)p^ + (l-0)p^_2 







[< 


(1-v^) 


Pi + 


i-1 


(27c) 


+ V, 


(l+0)p^_j^ + 


(1-0)P 


L-2 


(4>i-<^i-i) 


2 
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(27d) 


Here, A5 “ An “ 1 and only the varying indices are indicated. The param- 
eter 0 = 1 or 2 for first- or second-order spatial accuracy in super- 
sonic regions. The switching parameter v is defined in a way similar to 
(17) as 


V = ” (p/p*)^3c 1 ^ c < 10 

V = 0 if V < 0 (i.e., subsonic) 

V = 1 if V > 1 (i.e., supersonic) 


(28) 


where p* is density evaluated at sonic conditions. 

The parameter v can be set to 1 throughout, but accuracy will be 
Impaired unless 6 is also set to 2. The operators in Eqs. (27a) and 
(27b) assume that the flow will be supersonic only in the positive 
x-direction. The density is found from the Bernoulli equation with 
(AC = An = 1): 


. * 1+1 ■ * 1-1 




2 n j 


At 


The metrics 5 and n are obtained from 
X z 


6 ()> 
X 


n+1 


«x “ - x^_, 


■ 'j-1 

while the term (^6) is formed either as 






i+% 


(29a) 
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or 







(29b) 


The terms (5x/J)i_i/2> (nz/J)j+l/2» and (nz^/J)j-l /2 receive similar 
treatment. If Eq. (19a) is used, it is essential to add -6c(Pa,^ ^/j)Se$„ 
to the right-hand side of Eq. (16) to subtract out a numerical truncation 
error due to incomplete metric cancellation. 

Eq. 26 is now rearranged into delta form; it is now to be solved for 
- <|,n+l ^ ^n. For example, the term 6^(C^2 3n)6^^n+l can be rearranged 


When the equation is put into delta form all the unknown terms are 
arranged together, resulting in 

- h6^(n//J)“*‘ - ♦“) 

= + (■>/)" C’J,}*"-*"-') ( 

+ (e” - 0"‘S + h(s^(c2/J)'«-* p“6^+'’- 

which can be approximately factored into the form 
(I + at(„ 2)“+‘ ♦ “j - it(j-‘^V6")hs (n 2/J)”** p"5 )x 

7 n n M y f) 

{I + At(C^2)n+l ^n^^ _ At(j"+V8“)h6^(C^2/J)“+^p^6^}(^“+l-^“) 

= [1 + (6‘''*V6")(J"'^Vj")]( 4.“- - (6""Ve")(j“'^Vj“)(<|.“‘^ - 

+ At(s"'Vs")(j"‘^Vj") + (ny2)"4>"'^5j(<f“-<f."'b ^ 

+ At(j"+Ve"){(p" - p"’^) + h6.(52/J)"^^p"6 ()." 

^ X 5 

+ h6^(n2/j)”‘*‘^p^6^4,"} 
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This equation has the form 


- 4,") = R 

and is implemented as an algorithm as 


(28) 


L A(|)* = R 

n 






A({.* 

4,“ + A4i" 


(29) 


This algorithm requires only a series of scalar tridiagonal inversions 
and it is therefore very efficiently implemented. Computer storage equiv- 
alent to four levels of have to be supplied with p computed from the 
Bernoulli equation as needed. 

RESULTS AND APPLICATIONS 


Although the various finite difference issues are very similar for 
the unsteady small disturbance and full potential methods, there is a 
large difference in the extent to which these methods have been developed 
and applied. Naturally, this difference reflects the relative complexity 
of the two methods. 

One of the most Interesting test computations demonstrating the use 
of unsteady potential finite difference methods concerns the flow about 
a pulsating airfoil (that is, one having a time-varying thickness). 
Although the idea of a pulsating airfoil seems farfetched, it can be 
shown by considering the two-dimensional small disturbance equation to 
emulate a time-varying free-stream Mach number. And free-stream Mach 
number-variation does occur for an advancing helicopter rotor. Now con- 
sider a parabolic arc airfoil whose mid-chord thickness varies as 


/'0.1[ 10 - t + 6(t/15)2](t/15)3 , 0£t^l5 

T(« - { 0.l[l0 - is(^) + . 15 < t < 30 

0 , t > 30 


where t = t'Ueo/i and is nondlmensionallzed by chord lengths traveled. 
Since the variation takes place over many chords of travel it is suitable 
to invoke the low frequency approximation for the small disturbance equa- 
tion — that is, all time derivatives except are Ignored in Eq. (7). 

(Physically, ignoring in Eq. (7) amounts to assuming an infinite 

downstream propagation rate.) Figures 2 and 3 show a comparison of the 
resulting flow computed by the full potential method (Goorjian, Ref. 9) 
and the low frequency small disturbance equation (Ballhaus and Steger, 
Ref. 23). When the airfoil is thinning it becomes subcrltlcal by propa- 
gating the shock upstream from the leading edge. The two approaches 
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give essentially the same result. The most significant difference between 
the two computations is perhaps the greater dissipation of the forward 
propagating shock for the full potential case. This can be controlled by 
the choice of upstream density biasing function. For this case, density 
switching was employed. Recall that it is possible to use an unswltched 
upstream density as long as a higher order difference is used to maintain 
accuracy. This is demonstrated in Fig. 4 which compares switched and 
unswitched density biasing for the computation of the steady flow on a 
biconvex profile. 

The previously mentioned shock motion is so different from our 
previous experience that it demands some kind of experimental study. By 
coincidence, such a study was performed by Tlj deman and his associates 
(24,25) at NLR, the Netherlands, at the same time that the first of these 
computations was being done. They acquired detailed flow visualization 
and loading data for a NACA64A006 airfoil with an oscillating flap. They 
delineated three basic types of shock motion caused by the oscillating 
flap. These are: 

1. Type A. The shock moves nearly sinusoidally (only the lowest 
harmonic was measured) with a phase shift relative to the flap 
motion. The shock strength varies, being a minimum while moving 
downstream and a maximum while moving upstream. 

2. Type B. This case is similar to Type A except that the shock 
strength variation disappears during the downstream moving 
portion of its cycle. 

3. Type C. At slightly supercritical conditions there is no super- 
sonic region for a large portion of the flap cycle. In this 
case the airfoil becomes subcrltlcal by propagating the shock 
upstream off of the leading edge. There is no downstream shock 
motion. 

The ability of potential methods to compute these shock motions was 
demonstrated by Ballhaus and Goorjian (26) . In this work, two-dimensional 
finite difference solutions of the low frequency small disturbance equa- 
tion were obtained by the ADI approach (program, LTRAN2) . The computed 
Type B motion is shown in Fig. 5. The shock-motion disappearance and 
reappearance are shown here for LTRAN2 and the Magnus-Yoshihara code (an 
isentropic Euler code). Type C motion is illustrated in Fig. 6. In the 
succession of plots shown, the shock is seen to form at mid— chord and 
move upstream until it disappears at the leading edge. The shock is seen 
to disappear here rather than visibly propagate upstream as in Fig. 3, 
probably due to numerical dissipation. It is interesting that the condi- 
tions for which these various types of shock motion were computed do not 
compare well with the actual experimental conditions. It is generally 
felt now (on the basis of comparing different codes and of computations 
of wall effects) that the discrepancy is due to tunnel-wall and possibly 
viscous effects rather than numerical problems. Nevertheless, the fact 
that experiment and computation produce the same kinds of phenomena 
greatly increases the significance of both. 

A primary application of unsteady transonic potential methods is in 
the prediction of loads on helicopter rotors in forward flight. Although 
aeroelastic effects are important, in this case the main source of 
unsteadiness is in the flow itself. The most notable distinction between 
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the fixed and rotary wing is that for the latter the free stream (in body- 
fixed coordinates) is constantly accelerating and decelerating. That is, 
the rotor "sees" the free-stream Mach number to be periodically varying. 

The reduced frequency of this variation is the inverse of the blade 
aspect ratio. Since aspect ratio is of the order of 10 for a rotor, it is 
possible to invoke the low frequency approximation. The effectiveness of 
the small disturbance potential methods has been demonstrated by comparison 
of computed flows with measured surface pressure on a nonlifting rotor 
blade. Figures 7 and 8 demonstrate this comparison at two blade azimuths, 

= 60° and 120° , respectively (see Fig. 7 for definition of azimuth 
angle). Measurements and computations are in excellent agreement. The 
effects of unsteadiness are evident in this figure because for the two 
azimuths shown the chordwise Mach numbers are identical. Yet the pres- 
sures are quite different. There are no shocks seen at = 60“ , but they 
are evident at 120° . This azimuthal shock asymmetry is not explainable 
by cross-flow effects, because the inboard station is too far from the 
tip. Furthermore, the inboard results can be readily obtained by two- 
dimensional computations (27) . Very often, the shock motion in these 
rotor computations appears to be Type C, the upstream propagating type. 

This is seen in Fig. 9 (27) which shows a low frequency small disturbance 
two-dimensional computation of a lifting rotor flow. In this case the 
blade is oscillating and also sees a varying free-stream Mach number. 

The flow on the bottom surface of the airfoil is seen to return to sub- 
critical conditions by propagating the shock upstream. 

The most fascinating feature of the above computations is that they 
are so easy to perform compared with wind-tunnel testing. Yet they have 
often proven to contain most of the essential physics. In the following 
discussion we shall demonstrate the use of a potential computation as a 
"numerical wind tunnel" to explore a little-studied but possibly very 
important problem. 

An unusual unsteady flow feature of helicopter rotors is that they 
are never very far from the tip vortex of a preceding blade and close 
blade/vortex interactions often occur. Under certain conditions a 
blade can encounter a vortex which is nearly parallel to itself. Such 
encounters are an important noise source and require modeling. An 
extremely simple model is provided by a two-dimensional small distur- 
bance computation of a near-vortex encounter. The time scale for this 
problem (i/Uoo) is very brief and it is necessary to include all time 
derivatives. The vortex is introduced as the edge of a potential dis- 
continuity sheet (Fig. 10) which is stepped through the computational 
grid. The vortex is moved through the grid in a prescribed straight line 
at the undisturbed flow speed. The strength of the vortex is given as 
an effective lift coefficient, Clv» of an airfoil having the same circu- 
lation as the vortex. Note that while the sheet describing the vortex 
in Fig. 10 is horizontal, its direction is irrelevant. The effect of 
unsteadiness in these computations is shown in Fig. 11, which compares 
blade surface pressure distributions for a fixed and moving vortex (the 
vortex, located at mid-chord and 0.96 chords below the blade, has the 
strength, Clv ® 0.1). There is no apparent disturbance for the unsteady 
moving vortex case but a large disturbance for the steady case. This 
result undoubtedly reflects the fact that the vortex exerts no force on 
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the fluid In the unsteady case. In fact, in order to get a sizeable 
effect on the surface pressures (other than the expected angle-of.-attack 
variation) it is necessary to bring the vortex quite close to the blade 
and increase its strength. Figure 12 shows the lift variation for a 
blade with a vortex of strength Clv = 0.4 whose path lies 0.25 chords 
beneath the blade. As the vortex passes it induces a shock on the bottom 
surface. This shock disappears very suddenly with no obvious upstream 
propagation. Through all this activity the upper surface seems curiously 
unaffected. 

The importance of such a computational experiment is difficult to 
assess because many liberties have been taken. Nevertheless, it is clear 
that the neglect of unsteadiness would be to greatly overestimate the 
effects of these phenomena. For sufficiently strong or close vortices, 
large and rapid bottom surface disturbances can be computed and these 
must have acoustic significance. However, this could be mitigated by 
allowing the vortex to move freely rather than follow a fixed path. 
Undoubtedly, an assessment of these sorts of computations cannot be made 
without some experiments. Thus it seems that the effect of the "numerical 
wind tunnel" can be to guide the use of and increase the need for the 
physical wind tunnel. 

CONCLUDING REMARKS 

This paper has made no attempt to treat unsteady potential finite 
difference methods exhaustively or in great detail. Rather, a typical 
small disturbance method, and a full potential method have been discussed 
in parallel in the context of the issues (choice of equations and boundary 
conditions, linearizations, discretizations, and algorithms) which arise 
in all code development. One can discuss a typical method because the 
choices taken in code development are quite few. For instance, the choice 
i^^ the equations centers around whether one wishes to have one or two 
dependent variables. Linearizations and discretizations vary little. 

All methods employ some sort of upstream biasing and nearly the same time 
discretizations. Algorithm development always comes down to some sort of 
approximate factorization. In short, although there is much development 
work to be done, the general area of potential finite difference methods 
is beginning to mature quite nicely. We shall surely see the development 
of several practical three-dimensional unsteady full potential codes 
within the coming year or two. 

There remains one vital issue which has not been covered in this dis- 
cussion: the subject of grid generation. On this point finite difference 

methods are often more of an art than a science. The problem is not as 
much to generate a mesh as to know the effect of this mesh on a particular 
problem and solution method. This is especially difficult in the unsteady 
case where we may have the presence of moving flow features (shocks and 
vortices) which require resolution. It is ironic that while gridding is 
the most fundamental feature of finite difference methods, the topic of 
grids remains to be organized into an organic and systematic entity. The 
requirements for large-scale unsteady computations must certainly change 
this situation, because they will require fast, automatic, and reliable 
grid schemes. 
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Fig. 5 - Unsteady upper surface pressure variation for a 
NACA64A006 airfoil with an oscillating flap. Type B motion 
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Fig. 6 - Unsteady upper surface pressure variation for a 
NACA64A006 airfoil with an oscillating flap, Type C motion 
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Fig. 7 - Measured and computed pressure distribution on 
an advancing, nonlifting rotor — 60“ blade azimuth 
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Fig. 8 - Measured and computed pressure distribution on 
an advancing, nonlifting rotor — 120® blade azimuth 
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